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I . INTRODUCTION 

The Penn State Finite Difference Time Domain Electromagnetic 
Scattering Code Version D is a three dimensional numerical 
electromagnetic scattering code based upon the Finite Difference 
Time Domain Technique (FDTD) . The supplied version of the code 
is one version of our current three dimensional FDTD code set. 
This manual provides a description of the code and corresponding 
results for several scattering problems. The manual is organized 
into fourteen sections: introduction, description of the FDTD 

method, operation, resource requirements, Version D code 
capabilities, a brief description of the default scattering 
geometry, a brief description of each subroutine, a description 
of the include file (COMMOND. FOR) , a section briefly discussing 
Radar Cross Section (RCS) computations, a section discussing some 
scattering results, a sample problem setup section, a new problem 
checklist, references and figure titles. 

II. FDTD METHOD 

The Finite Difference Time Domain (FDTD) technique models 
transient electromagnetic scattering and interactions with 
objects of arbitrary shape and/or material composition. The 
technique was first proposed by Yee [1] for isotropic, non- 
dispersive materials in 1966; and has matured within the past 
twenty years into a robust and efficient computational method. 

The present FDTD technique is capabable of transient 
electromagnetic interactions with objects of arbitrary and 
complicated geometrical shape and material composition over a 
large band of frequencies. This technique has recently been 
extended to include dispersive dielectric materials, chiral 
materials and plasmas. 

In the FDTD method, Maxwell's curl equations are discretized 
in time and space and all derivatives (temporal and spatial) are 
approximated by central differences. The electric and magnetic 
fields are interleaved in space and time and are updated in a 
second-order accurate leapfrog scheme. The computational space 
is divided into cells with the electric fields located on the 
edges and the magnetic fields on the faces (see Figure 1) . FDTD 
objects are defined by specifying dielectric and/or magnetic 
material parameters at electric and/or magnetic field locations. 

Two basic implementations of the FDTD method are widely used 
for electromagnetic analysis: total field formalism and scattered 
field formalism. In the total field formalism, the electric and 
magnetic field are updated based upon the material type present 
at each spatial location. In the scattered field formalism, the 
incident waveform is defined analytically and the scattered field 
is coupled to the incident field through the different material 
types. For the incident field, any waveform, angle of incidence 
and polarization is possible. The separation of the incident and 



scattered fields conveniently allows an absorbing boundary to be 
employed at the extremities of the discretized problem space to 
absorb the scattered fields. 
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This code is a scattered field code, and the total E and H 
fields may be found by combining the incident and scattered 
fields. Any type of field quantity (incident, scattered, or 
total) , Poynting vector or current are available anywhere within 
the computational space. These fields, incident, scattered and 
total, may be found within, on or about the interaction object 
placed in the problem space. By using a near to far field 
transformation, far fields can be determined from the near fields 
within the problem space thereby affording radiation patterns and 
RCS values. The accuracy of these calculations is typically 
within a dB of analytic solutions for dielectric and magnetic 
sphere scattering. Further improvements are expected as better 
absorbing boundary conditions are developed and incorporated. 

III. OPERATION 


Typically, a truncated Gaussian incident waveform is used to 
excite the system being modeled, however certain code versions 
also provide a smooth cosine waveform for convenience in modeling 
dispersive materials. The interaction object is defined in the 
discretized problem space with arrays at each cell location 
created by the discretization. All three dielectric material 
types for E field components within a cell can be individually 
specified by the arrays IDONE (I , J, K) , IDTWO(I, J,K) , 

IDTHRE (I , J, K) . This models arbitrary dielectric materials with /x 
= Mo* B Y an obvious extension to six arrays, magnetic materials 
with /x f n 0 can be modeled. 


Scattering occurs when the incident wave, marched forward in 
time in small steps set by the Courant stability condition, 
reaches the interaction object. Here a scattered wave must 
appear along with the incident wave so that the Maxwell equations 
are satisfied. If the material is a perfectly conductive metal 
then only the well known boundary condition 


scat inc 

Etan = "B tan 


(1) 


must be satisfied. For a nondispersive dielectric the 
requirement is that the total field must satisfy the Maxwell 
equations in the material: 


Vx E tot 


Vx (E inc + E scat 


1 8H tot 
M 0 dt 


1 d(H inc +H scat ) 


at 


(2) 
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a-ptot 

Vx H tot = Vx (H inc +H scat ) =e — — — +aE tot 

at 


(3) 


= e 


a(E ,nc +E ! 


scat 


at 


- +a (E inc + E scat ) 


(4) 


Additionally the incident wave, defined as moving unimpeded 
through a vacuum in the problem space, satisfies everywhere in 
the problem the Maxwell equations for free space 


Vx E inc 


i aH inc 
M 0 dt 


(5) 


VxH inc 



( 6 ) 


Subtracting the second set of equations from the first yields the 
Maxwell equations governing the scattered fields in the material: 


Vx E scat 


i aH scat 
at 


(7) 


ap' nc ap scat 

Vx H scat = (e -e ) — — +aE ,nc +e — + crE scat 

0 at at 


( 8 ) 


Outside the material this simplifies to: 

Vx E scat = - — 


i aH scat 


M 0 at 


(9) 


VxH scat =e 


dE 


scat 


at 


( 10 ) 


Magnetic materials, dispersive effects, non-linearities, 
etc. , are further generalizations of the above approach. Based 
on the value of the material type, the subroutines for 
calculating scattered E and H field components branch to the 
appropriate expression for that scattered field component and 
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that component is advanced in time according to the selected 
algorithm. As many materials can be modeled as desired, the 
number equals the dimension selected for the flags. If materials 
with behavior different from those described above must be 
modeled, then after the appropriate algorithm is found, the 
code's branching structure allows easy incorporation of the new 
behavior. 

IV. RESOURCE REQUIREMENTS 

The number of cells the problem space is divided into times 
the six components per cell set the problem space storage 
requirements 

Storage = NC x 6 components/cell x 4 bytes/component (11) 
and the computational cost 

Operations = NC x 6 comp/cell x 10 ops/component x n (12) 

where N is the number of time steps desired. 

N typically is on the order of ten times the number of cells 
on one side of the problem space. More precisely for cubical 

cells it takes v^3 time steps to traverse a single cell when the 
time step is set by the Courant stability condition 

, Ax . 

At = Ax = cell size dimension 

/7c 


The condition on N is then that 

1 1 

N - lOx (/7 nC^) NC* 


number cells on a side 
of the problem space 


(14) 


The earliest aircraft modeling using FDTD with approximately 30 
cells on a side required approximately 500 time steps. For more 
recent modeling with approximately 100 cells on a side, 2000 or 
more time steps are used. 

For (100 cell) 3 problem spaces, 24 MBytes of memory are 
required to store the fields. Problems on the order of this size 
have been run on a Silicon Graphics 4D 220 with 32 MBytes of 
memory, IBM RISC 6000, an Intel 486 based machine, and VAX 
11/785. Storage is only a problem as in the case of the 486 
where only 16 MBytes of memory was available. 3 This limited the 
problem space size to approximately (80 cells) . 
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For (100 cell) problems with approximately 2000 time steps, 
there is a total of 120 x 10 7 8 9 operations to perform. The speeds 
of the previously mentioned machines are 24 MFLOPs (4 processor 
upgraded version), 10 MFLOPS, 1.5 MFLOPS, ^nd 0.2 MFLOPs. T^e 
run times are then 5 3 x 10 3 seconds, 12 x 10 3 seconds, 80 x 10 3 
seconds and 600 x 10 3 seconds, respectively. In hours the times 
are 1.4, 3.3, 22.2 and 167 hours. Problems of this size are 
possible on all but the last machine and can in fact be performed 
on a personal computer (486) if one day turnarounds are 
permissible. 

V. VERSION D CODE CAPABILITIES 

The Penn State University FDTD Electromagnetic Scattering 
Code Version D has the following capabilities: 

1) Ability to model lossy dielectric and perfectly conducting 
scatterers . 

2) Ability to model lossy magnetic scatterers. 

3) Ability to model dispersive dielectric and dispersive 
magnetic scatterers. This dispersive FDTD method is now 
designated (FD) TD for Frequency-Dependent Finite Difference Time 
Domain. 

4) First and second order outer radiation boundary condition 
( ORBC ) operating on the electric fields for dielectric and 
dispersive dielectric scatterers. 

5) First and second order ORBC operating on the magnetic fields 
for magnetic and dispersive magnetic scatterers. 

6) Near to far zone transformation capability to obtain far zone 
scattered fields. 

7) Gaussian and smooth cosine incident waveforms with arbitrary 
incidence angles. 

8) Near zone field, current or power sampling capability. 

9) Companion code for computing Radar Cross Section (RCS) . 

VI. DEFAULT SCATTERING GEOMETRY 

The code as delivered is set up to calculate the far zone 
backscatter fields for a 6.67 meter radius, dispersive, Nickel 
Ferrite sphere. Nickel Ferrite is defined by a frequency 
dependent permeability given by 

where is the infinite frequency permeability, is the zero 



9 


M ( (ji) M s M ffl 

= fj, + 

M 0 " l + j«7 0 


(15) 


frequency permeability, r Q is the relaxation time and o> is the 
radian frequency. The Nickel Ferrite parameters are n n =l, n s =100 
and t q =20 ns. The problem space size is 66 by 66 by 66 cells in 

the x, y and z directions, the cells are 1/3 m cubes, and the 
incident waveform is a 0-polarized smooth cosine pulse with 
incidence angles of 0=22.5 and 0=22.5 degrees. The output data 
files are included as a reference along with a code (RCS3D. FOR) 
for computing the frequency domain RCS using these output data 
files. The ORBC is the second order absorbing boundary condition 
set forth by Mur [2]. 

VII. SUBROUTINE DESCRIPTION 

In the description for each subroutine, an asterisk (*) will 
be placed by the subroutine name if that particular subroutine is 
normally modified when defining a scattering problem. 

MAIN ROUTINE 

The main routine in the program contains the calls for all 
necessary subroutines to initialize the problem space and 
scattering object (s) and for the incident waveform, far zone 
transformation, field update subroutines, outer radiation 
boundary conditions and field sampling. 

The main routine begins with the include statement and then 
appropriate data files are opened, and subroutines ZERO, BUILD 
and SETUP are called to initialize variables and/or arrays, build 
the object (s) and initialize the incident waveform and 
miscellaneous parameters, respectively. Subroutine SETFZ is 
called to intialize parameters for the near to far zone 
transformation if far zone fields are desired. 

The main loop is entered next, where all of the primary 
field computations and data saving takes place. During each time 
step cycle, the EXSFLD , EYSFLD , and EZSFLD subroutines are called 
to update the x, y, and z components of the scattered electric 
field. The six electric field outer radiation boundary 
conditions (RADE??) are called next to absorb any outgoing 
scattered fields for perfectly conducting, dielectric, or 
dispersive dielectric scatterers. Time is then advanced 1/2 time 
step according to the Yee algorithm and then the HXSFLD, HYSFLD, 
AND HZSFLD subroutines are called to update the x, y, and z 
components of scattered magnetic field. The six magnetic field 
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outer radiation boundary conditions (RADH??) are called next to 
absorb any outgoing scattered fields for magnetic or dispersive 
magnetic scatterers. Time is then advanced another 1/2 step and 
then either near zone fields are sampled and written to disk in 
DATSAV, and/or the near zone to far zone vector potentials are 
updated in SAVFZ. The parameter NZFZ (described later) in the 
common file defines the type of output fields desired. 

After execution of all time steps in the main field update 
loop, subroutine FAROUT is called if far zone fields are desired 
to compute the far zone fields and write them to disk. At this 
point, the execution is complete. 

SUBROUTINE SETFZ 

This subroutine initializes the necessary parameters 
required for far zone field computations. The code as furnished 
computes backscatter far zone fields and can compute bistatic far 
zone fields for one scattering angle (i.e. one 0 and <p angle). 
Refer to reference [3] for a complete description of the near to 
far zone transformation. Other versions of this subroutine 
provide for multiple bistatic angles. 

SUBROUTINE SAVFZ 

This subroutine updates the near zone to far zone vector 
potentials. 

SUBROUTINE FAROUT 

This subroutine changes the near zone to far zone vector 
potentials to far zone electric field 0 and <p components and 
writes them to disk. 

SUBROUTINE BUILD * 

This subroutine "builds” the scattering object (s) by 
initializing the I DONE, IDTWO, IDTHRE, IDFOR, IDFIV and IDSIX 
arrays. The I DONE - I DTHRE arrays are for specifying perfectly 
conducting, lossy dielectrics and dispersive dielectric 
materials. The IDFOR-IDSIX arrays are for lossy magnetic and 
dispersive magnetic materials. The reason for the separate 
arrays is so the user can independently control the exact 
placement of dielectric and magnetic material in the Yee cells. 
Refer to Figure 1 for a diagram of the basic Yee cell. For 
example, setting an element of the I DONE array at some I,J,K 
location is actually locating dielectric material at a cell edge 
whose center location is I+0.5,J,K. Setting an element of the 
IDFOR array at some I,J,K location is actually locating magnetic 
material perpendicular to a cell face whose center location is 
I, J+0. 5,K+0. 5, or equivalently, locating magnetic material at an 
edge on the dual H field mesh. The difference between the IDONE 
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and IDFOR array locations is a direct result of the field offsets 
in the Yee cell (see Figure 1) . Thus, materials with diagonal 
permittivity and/or diagonal permeability tensors can be modeled. 
The default material type for all ID??? arrays is 0, or free 
space. By initializing these arrays to values other than 0, the 
user is defining an object by determining what material types are 
present at each spatial location. Other material types available 
for IDONE-IDTHRE are 1 for perfectly conducting objects, 2-9 for 
lossy non-magnetic dielectrics, 20-29 for dispersive dielectrics. 
IDONE-IDTHRE are normally set to 0 for magnetic scatterers. Other 
material types available for IDFOR-IDSIX are 10-19 for lossy 
magnetic materials and 30-39 for dispersive magnetic materials. 
IDFOR-IDSIX are normally set to 0 for perfectly conducting or 
dielectric scatterers. If the user wants a material with both 
dielectric and magnetic properties (i.e. permittivity other than 
e 0 for magnetic materials, and permeability other than n 0 for 
dielectric materials) , then he/ she must define IDONE-IDSIX for 
that particular material. This subroutine also has a section 
that checks the ID??? arrays to determine if legal material types 
have been defined throughout the problem space^. The actual non- 
dispersive material parameters (e, n, a, and a ) are defined in 

subroutine SETUP. The dispersive material parameters (e s , 

t 0 , a, / t q , and a*) are also defined in a separate section 

in SETUP. The default geometry is a 6.67 m radius, dispersive, 
Nickel Ferrite sphere. 

The user must be careful that his/her object created in the 
BUILD subroutine is properly formed. There is not a direct 
one-to-one correspondence between the dielectric and magnetic 
ID??? arrays. However, one can define a correspondence, so that 
code used to generate a dielectric object can be modified to 
generate a magnetic object. 

To see this consider that we have set the permittivity at 
cell locations corresponding to 

EX ( I , J , K) , EY ( I , J , K) , EZ ( I , J, K) 

using the IDONE, IDTWO, and IDTHRE arrays respectively. This 
determines one corner of a dielectric cube. If we wish to define 
the corner of a corresponding magnetic cube, offset 1/2 cell in 
the x, y, z directions, we would set the locations of the 
magnetic fields 

HX ( 1+1 , J , K) , HY ( I , J+l , K) , HZ ( I , J, K+l ) 
as magnetic material using the IDFOR, IDFIV, and IDSIX arrays. 

This example indicates the following correspondence between 
BUILDing dielectric and magnetic objects: 
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Dielectric Magnetic 

Object Object 

IDONE ( I , J , K) = IDFOR ( 1+1 , J , K) 

IDTWO ( I , J , K) = IDFIV ( I , J+l , K) 

IDTHRE (I , J, K) = IDSIX ( I , J , K+l ) 

To illustrate this for a somewhat more complicated case, consider 
an example of a 3x3 cell dielectric plate located in the XY plane 
at K=5. The plate can be generated with the following FORTRAN 
lines: 

PLATE 

7. . 


6 

J 

A 5 


4 . • 

4 5 6 7 

>x 

If the same FORTRAN lines were used to try to generate a magnetic 
plate, the object generated would actually be unconnected: 

PLATE 

7. . 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF (I • NE. 7 ) IDONE (I, J,K)=2 
IF ( J .NE. 7) IDTWO (I, J,K)=2 
10 CONTINUE 
20 CONTINUE 


PLATE 


i 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF (I .NE. 7) IDFOR (I , J, K) =11 
IF ( J . NE . 7 ) IDFIV ( I , J , K) =11 
10 CONTINUE 
20 CONTINUE 




The correct way to build the magnetic plate is with the following 
FORTRAN code : 


DO 20 J=4 , 7 
DO 10 1=4,7 

IF ( I . NE . 4 ) IDFOR (I,J,K)=11 

IF ( J . NE . 4 ) IDFIV (I , J, K) =11 
10 CONTINUE 
20 CONTINUE 
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This is equivalent to the correspondence relationship given 
above. See comments in the BUILD subroutine for further 
explanation of the ID??? arrays. 

When it is important to place the object in the center of 
the problem space (to have lowest possible cross-pol scattering 
for symmetric objects) , NX etc. should be odd for dielectrics and 
even for magnetics. This is due to the field locations in the 
Yee cell and also the placement of the E and H field absorbing 
boundary condition surfaces. 

If the object being modeled has curved surfaces, edges, etc. 
that are at an angle to one or more of the coordinate axes, then 
that shape must be approximately modeled by lines and faces in a 
"stair-stepped” (or stair-cased) fashion. This stair-cased 
approximation introduces errors into computations at higher 
frequencies. Intuitively, the error becomes smaller as more 
cells are used to stair-case a particular object. Thus, the 
default Nickel Ferrite sphere scattering geometry is a stair- 
cased sphere. 

When the user's test object is dielectric it is apparently 
better to use the Mur E field boundary formulation. For a 
magnetic scattering object it seems better to use the Mur H field 
boundary formulation. The reason for this has to do with the 
reactive part of the energy (near zone fields) reaching the 
absorbing boundary surface. These near zone fields differ for 
magnetic and dielectric objects, even for dual objects where the 
far fields are essentially identical. 

One final comment on building dual dielectric vs magnetic 
objects can be explained by considering the duality between 
magnetic and dielectric materials. In testing the code, it was 
helpful to test dual dielectric/magnetic objects since they 
scatter identically in the far zone. In specifying the material 
of a magnetic scatterer to be the dual of a dielectric, the 
following transformation must be applied: 


DIELECTRIC 

< — DUAL- 

-> 

MAGNETIC 

E field 

< 

-> 

H field 

H field 

< 

-> 

-E field 

e 

< 

-> 

M 

M 

< 

-> 

£ 

fl (k) 

< 

-> 

A ( k > 

a 

< 

-> 

* (M 0 /e 0 ) =a 


The somewhat surprising entry in this table is the relationship 
between the dielectric and magnetic conductivities. The reason 
why the magnetic conductivity, a , for magnetic materials has to 
be scaled as indicated is that for duality to be applied the free 
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space impedance must be inverted. Although this is not generally 
possible for actual problems, identical far zone scattering for 
dielectric apd magnetic scatterers can be achieved by this 
scaling of a as above and realizing that E and H scattered £ield 
to incident field ratios will not invert. This scaling of a is 
why conductivity for magnetic materials is usually not a 
dominating feature, and in fact is often neglected. 

SUBROUTINE DCUBE 

This subroutine builds cubes of dielectric material by 
defining four each of I DONE, IDTWO and IDTHRE components 
corresponding to one spatial cube of dielectric material. It can 
also be used to define thin (i.e. up to one cell thick) 
dielectric or perfectly conducting plates. Refer to comments 
within DCUBE for a description of the arguments and usage of the 
subroutine. This subroutine could be modified to build cubes 
and/or plates of magnetic materials by using a triple do-loop in 
BUILD (after calls to DCUBE) over coordinate indices I, J, K and 
applying the correspondence between the IDONE-IDTHRE and IDFOR- 
IDSIX arrays that was previously discussed. 

SUBROUTINE SETUP * 

This subroutine initializes many of the constants required 
for incident field definition, field update equations, outer 
radiation boundary conditions and 4 material parameters. The 
material parameters e, fi, a and a are defined for each material 
type (non-dispersive) using the material arrays EPS, XMU, SIGMA 
and SIGMAC respectively. The array EPS is used for the total 
permittivity, XMU is used for the total permeability, SIGMA is 
used for the electric conductivity and SIGMAC is used for the 
magnetic conductivity (useful for running dual problems) . These 
arrays are initialized in SETUP to free space material parameters 
for all material types and then the user is required to modify 
these arrays for his/her scattering materials. Thus, for the 
lossy dielectric material type 2, the user must define EPS(2) and 
SIGMA ( 2 ) . 

For dispersive dielectric and magnetic materials, different 
material parameter arrays are used. The functional form of the 
frequency dependent permittivity/permeability that was 
implemented in the code is the Debye relaxation [4] with an 
effective DC conductivity given by 

£(co) = e / -je // = e o e 0 + e 0 X e (<*>)+ (16) 

D<*> 


where the frequency dependent electric susceptibility function is 
defined as 
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X e (w) 


e g - e. 

1 + j'UX 0 


(17) 


where e a is the infinite frequency permittivity, e s is the zero 

frequency permittivity, t q is the relaxation time, a is the 

effective electric conductivity, and ca is the radian frequency. 
The same expressions were used for the frequency dependent 
permeability and can be written from equations (16) and (17) by 
replacing e by m and a by a . The corresponding time domain 
susceptibility function is given by 


x e (t) 



(18) 


The FDTD implementation of frequency dependent permittivity 
and/or permeability involves a convolution with the electric 
and/or magnetic field and interested readers are referred to 
references [5-6] for further details. 

For dispersive dielectric materials, the corresponding 
material parameter arrays are EPSINF (€„), EPSSTA (e s ), RELAXT 

(r Q ), and RELSIG (ct) . For dispersive magnetic materials, the 

arrays are XMUINF (/xj, XMUSTA (/x s ) , RELAXT ( t q ) and RELSIG (a*). 

These dispersive material parameters are defined under the 
DISPERSIVE SETUP portion of the subroutine. The remainder of the 
subroutine computes constants used in field update equations and 
boundary conditions and writes the diagnostics file. 

SUBROUTINE EXSFLD 

This subroutine updates all x components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDONE array are used to determine the type of material present 
and the corresponding update equation to be used. These 
scattered field equations are based on the development given in 
[7]. 

SUBROUTINE EYSFLD 

This subroutine updates all y components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTWO array are used to determine the type of material present 



and the corresponding update equation to be used. 
SUBROUTINE EZSFLD 
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This subroutine updates all z components of scattered 
electric field at each time step except those on the outer 
boundaries of the problem space. IF statements based upon the 
IDTHRE array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADEYX, RADEZX, RADEZY, RADEXY, RADEXZ and RADEYZ 

These subroutines apply the outer radiation boundary 
conditions to the scattered electric field on the outer 
boundaries of the problem space for non-magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 

SUBROUTINE HXSFLD 

This subroutine updates all x components of scattered 
magnetic field at each time step. IF statements based upon the 
IDFOR array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE HYSFLD 

This subroutine updates all y components of scattered 
magnetic field at each time step. IF statements based upon the 
IDFIV array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINE HZSFLD 

This subroutine updates all z components of scattered 
magnetic field at each time step. IF statements based upon the 
IDSIX array are used to determine the type of material present 
and the corresponding update equation to be used. 

SUBROUTINES RADHXZ, RADHYX, RADHZY , RADHXY , RADHYZ and RADHZX 

These subroutines apply the outer radiation boundary 
conditions to the scattered magnetic field on the outer 
boundaries of the problem space for magnetic scatterers. The 
user controls selection of these routines through the parameter 
MAGNET (described later) . 


SUBROUTINE DATSAV * 
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This subroutine samples near zone scattered field quantities 
and saves them to disk. This subroutine is where the quantities 
to be sampled and their spatial locations are to be specified and 
is only called if near zone fields only are desired or if both 
near and far zone fields are desired. Total field quantities can 
also be sampled. See comments within the subroutine for 
specifying sampled scattered and/or total field quantities. When 
sampling magnetic fields, remember the <5t/2 time difference 
between E and H when writing the fields to disk. Sections of 
code within this subroutine determine if the sampled quantities 
and the spatial locations have been properly defined. 

FUNCTIONS EXI, EYI , EZI, HXI , HYI and HZI 

These functions are called to compute the x, y and z 
components of incident electric and magnetic field, respectively. 
The functional form of the incident field is contained in a 
separate function SOURCE. 

FUNCTION SOURCE * 

This function contains the functional form of the incident 
field. The code as furnished uses the smooth cosine form of the 
incident field. An incident Gaussian pulse is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, this function need only be modified if the 
user changes the incident pulse from smooth cosine to Gaussian. 
Currently, only the smooth cosine pulse can be used for 
scattering from dispersive targets. A slight improvement in 
computing speed and vectorization may be achieved by moving this 
function inside each of the incident field functions EXI, EYI and 
so on. 

FUNCTIONS DEXI, DEYI , DEZI, DHXI , DHYI and DHZI 

These functions are called to compute the x, y and z 
components of the time derivative of incident electric and 
magnetic field, respectively. The functional form of the 
incident field is contained in a separate function DSRCE. 

FUNCTIONS DEXIXE, DEYIXE, DEZIXE, DHXIXE, DHYIXE and DHZIXE 

These functions compute the x, y and z components of the 
convolution of the time derivative of the incident field with the 

electric or magnetic Debye susceptibility function x • 


FUNCTION DSRCE * 
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This function contains the functional form of the time 
derivative of the incident field. The code as furnished uses the 
time derivative of the smooth cosine form of the incident field. 

A Gaussian pulse time derivative is also available by 
uncommenting the required lines and commenting out the smooth 
cosine pulse. Thus, the function need only be modified if the 
user changes from the smooth cosine to Gaussian pulse. Again, a 
slight improvement in computing speed and vectorization may be 
achieved by moving this function inside each of the time 
derivative incident field functions DEXI, DEYI and so on. 

FUNCTION DCONV 

This function evaluates the convolution of the time 
derivative of the incident field with the Debye susceptibility 

function x • 

SUBROUTINE ZERO 

This subroutine initializes various arrays and variables to 

zero. 

VIII. INCLUDE FILE DESCRIPTION ( COMMOND . FOR ) * 

The include file, COMMOND. FOR, contains all of the arrays 
and variables that are shared among the different subroutines. 
This file will reguire the most modifications when defining 
scattering problems. A description of the parameters that are 
normally modified follows. 

The parameters NX, NY and NZ specify the size of the problem 
space in cells in the x, y and z directions respectively. For 
problems where it is crucial to center the object within the 
problem space, then NX, NY and NZ should be odd for dielectric 
scatterers and even for magnetic scatterers. The parameter NTEST 
defines the number of near zone guantities to be sampled and NZFZ 
defines the field output format. Set NZFZ=0 for near zone fields 
only, NZFZ=1 for far zone fields only and NZFZ=2 for both near 
and far zone fields. Parameter MAGNET is used to define magnetic 
scatterers and it controls the choice of RADE?? versus RADH?? 
absorbing boundary subroutines. It is set to 0 for dielectric 
scatterers and to 1 for magnetic scatterers. Parameter NUMMAT 
defines the total number of material types that are available for 
use. NEDISP and NHDISP define the number of dispersive 
dielectric and magnetic materials that are being used. Parameter 
NSTOP defines the maximum number of time steps. DELX, DELY, and 
DELZ (in meters) define the cell size in the x, y and z 
directions respectively. The 0 and 0 incidence angles (in 
degrees) are defined by THINC and PHINC respectively and the 
polarization is defined by ETHINC and EPHINC. ETHINC=1.0, 
EPHINC=0.0 for 6-polarized incident field and ETHINC=0.0, 
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EPHINC=1.0 for ^-polarized incident fields. Parameters AMP and 
BETA define the maximum amplitude and the e'^ temporal width of 
the incident pulse respectively. BETA automatically adjusts when 
the cell size is changed and normally should not be changed by 
the user. The far zone scattering angles are defined by THETFZ 
and PHIFZ. The code as furnished performs backscatter 
computations, but these parameters could be modified for a 
bistatic computation. 

IX. RCS COMPUTATIONS 

A companion code, RCS 3D. FOR, has been included to compute 
RCS versus frequency. It uses the file name of the (FD) TD far 
zone output data (FZ0UT3D.DAT) and writes data files of far zone 
electric fields versus time (FZTIME.DAT) and RCS versus frequency 
(3DRCS.DAT) . The RCS computations are performed up to the 10 
cell/A 0 frequency limit. Refer to comments within this code for 
further details. 


X. RESULTS 


As previously mentioned, the code as furnished models a 6.67 
m radius, dispersive. Nickel Ferrite sphere and computes 
backscatter far zone scattered fields at angles of 0=22.5 and 
0=22.5 degrees. Results are included for 0.25 dB loaded 
dielectric and magnetic foam spheres, 60 dB dielectric and 
magnetic foam spheres and the Nickel Ferrite magnetic sphere. 

The material parameters for 0.25 dB and 60 dB loaded foam are: 


0.25 DB FOAM 


60 DB FOAM 


1.16 e s 

1-01 e o 

0.6497 ns r Q 

2 . 954E-04 S/m <7 


41.0 

1.6 

0.3450 ns 
3 . 902E-01 S/m 


The magnetic materials are duals of the above as described 
earlier with the relative conductivity scaled by the ratio m o/ e o* 
The Nickel Ferrite parameters are ii s =100 , r Q =20 ns and 

a =0.0. For the 0.25 dB foams there are 10 cells per wavelength 
at approximately 3.0 GHz. For the 60 dB foam there are 10 cells 
per wavelength at approximately 2.35 GHz, as the 60 dB foam has a 
higher dielectric constant. 


Figures 2-10 show the real and imaginary parts of the 0.25 
dB and 60 dB foam permeability, the real and imaginary parts of 


the Nickel Ferrite permeability, and the magnetic 
susceptibilities versus time for all three materials. 
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Figures 11-14 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 0.25 dB and 60 dB 
dielectric foam spheres respectively. 

Figures 15-18 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the 0.25 dB and 60 dB 
mangetic foam spheres respectively. 

Figures 19-20 show the co-polarized far zone electric field 
versus time and the co-polarized RCS for the default Nickel 
Ferrite sphere using the dispersive FDTD method and the non- 
dispersive FDTD method. For the non-dispersive method, an 
equivalent permeability and magnetic conductivity were defined at 
3 0 MHz from (16) with /u replacing e as 

~ ^ I (o>=2t30E6) ' ° = ^ ^0 I (g>=2*30E6) 


XI. SAMPLE PROBLEM SETUP 

The code as furnished models a 6.67 m radius Nickel Ferrite 
dispersive magnetic sphere and computes backscatter far zone 
scattered fields at angles of 0=22.5 and 0=22.5 degrees. The 
corresponding output data files are also provided, along with a 
code to compute Radar Cross Section using these data files. In 
order to change the code to a new problem, many different 
parameters need to be modified. A sample problem setup will now 
be discussed. 

Suppose that the problem to be studied is RCS backscatter 
versus frequency from a 28 cm by 31 cm perfectly conducting plate 
with a 3 cm dielectric coating with a dielectric constant of 4e 0 
using a 0-polarized field. The backscatter angles are 0=30.0 and 
0=60 . 0 degrees and the frequency range is up to 3 Ghz . 

Since the frequency range is up to 3 Ghz, the cell size must 
be chosen appropriately to resolve the field IN ANY MATERIAL at 
the highest frequency of interest. A general rule is that the 
cell size should be 1/10 of the wavelength at the highest 
frequency of interest. For difficult geometries, 1/20 of a 
wavelength may be necessary. The free space wavelength at 3 GHz 
is A. 0 =lO cm and the wavelength in the dielectric coating at 3 GHz 
is 5 cm. The cell size is chosen as 1 cm, which provides a 
resolution of 5 cells/X in the dielectric coating and 10 cells/A, 0 
in free space. Numerical studies have shown that choosing the 

cell size < 1/4 of the shortest wavelength in any material is 
the practical lower limit. Thus the cell size of 1 cm is barely 
adequate. The cell size in the x, y and z directions is set in 
the common file through variables DELX, DELY and DELZ. Next the 
problem space size must be large enough to accomodate the 
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scattering object, plus at least a five cell boundary (10 cells 
is more appropriate) on every side of the object to allow for the 
far zone field integration surface. It is advisable for plate 
scattering to have the plate centered in the x and y directions 
of the problem space in order to reduce the cross-polarized 
backscatter and to position the plate low in the z direction to 
allow strong specular reflections multiple encounters with the 
ORBC . A 10 cell border is chosen, and the problem space size is 
chosen as 49 by 52 by 49 cells in the x, y and z directions 
respectively. As an initial estimate, allow 2048 time steps so 
that energy trapped within the dielectric layer will radiate. 

Thus parameters NX, NY and NZ in COMMOND.FOR would be changed to 
reflect the new problem space size, and parameter NSTOP is 
changed to 2048. If all transients have not been dissipated 
after 2048 time steps, then NSTOP will have to be increased. 
Truncating the time record before all transients have dissipated 
will corrupt frequency domain results. Parameter NZFZ must be 
equal to 1 since we are interested in far zone fields only. 
Parameter MAGNET must be equal to 0 for the dielectric scatterer. 
To build the object, the following lines are inserted into the 
BUILD subroutine: 

C 

C BUILD THE DIELECTRIC SLAB FIRST 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=3 1 

NZWIDE=3 

MTYPE=2 

CALL DCUBE ( I START , JSTART , KSTART , NXWIDE , NYWIDE , NZWIDE , MTYPE) 
C 

C BUILD PEC PLATE NEXT 

C 

ISTART=11 

JSTART=11 

KSTART=11 

NXWIDE=28 

NYWIDE=31 

NZWIDE=0 

MTYPE=1 

CALL DCUBE ( ISTART , JSTART , KSTART , NXWIDE , NYWIDE , NZWIDE , MTYPE ) 

The PEC plate is built last on the bottom of the dielectric 
slab to avoid any air gaps between the dielectric material and 
the PEC plate. In the common file, the incidence angles THINC 
and PHINC have to be changed to 30.0 and 60.0 respectively, the 
cell sizes ( DELX , DELY, DELZ) are set to 0.01, and the 
polarization is set to ETHINC=1.0 and EPHINC=0.0 for 0-polarized 
fields. Since dielectric material 2 is being used for the 
dielectric coating, the constitutive parameters EPS (2) and 



SIGMA ( 2 ) are set to 4e 0 and 0.0 respectively, in subroutine 
SETUP. This completes the code modifications for the sample 
problem. 

XII. MEW PROBLEM CHECKLIST 

This checklist provides a quick reference to determine if 
all parameters have been defined properly for a given scattering 
problem. A reminder when defining quantities within the code: 
use MKS units and specify all angles in degrees. 

COMMOND . FOR : 

1) Is the problem space sized correctly? (NX, NY, NZ) 

2) For near zone fields, is the number of sample points correct? 
(NTEST) 

3) Is parameter NZFZ defined correctly for desired field 
outputs? 

4) Is parameter MAGNET defined correctly for the type of 
scatterer? 

5) Is the number of dispersive dielectric (NEDISP) and 
dispersive magnetic (NHDISP) materials defined correctly? 

6) Is the number of time steps correct? (NSTOP) 

7) Are the cell dimensions (DELX, DELY, DELZ) defined correctly? 

8) Are the incidence angles (THINC, PHINC) defined correctly? 

9) Is the polarization of the incident wave defined correctly 
(ETHINC, EPHINC)? 

10) For other than backscatter far zone field computations, are 
the scattering angles set correctly? (THETFZ , PHIFZ) 

SUBROUTINE BUILD: 

1) Is the object completely and correctly specified? 

SUBROUTINE SETUP: 

1) Are the constitutive parameters for each material specified 
correctly? (EPS, XMU, SIGMA, SIGMAC) 


2) Are the constitutive parameters for each dispersive material 
defined correctly? (EPSSTA, EPSINF, RELAXT , RELSIG, XMUINF, 
XMUSTA, RELAXT, RELSIG) 
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FUNCTIONS SOURCE and DSRCE: 

1) If the smooth cosine pulse is not desired, is it commented 

out and the Gaussian pulse uncommented? 

SUBROUTINE DATSAV : 

1) For near zone fields, are the sampled field types and spatial 

locations correct for each sampling point? (NTYPE, IOBS, JOBS, 

KOBS) 
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XI. FIGURE TITLES 

Fig. 1 Standard three dimensional Yee cell showing placement 
of electric and magnetic fields. 

Fig. 2 Real part of relative permeability versus frequency for 
0.25 dB magnetic foam. 



24 


Fig. 3 
Fig. 4 
Fig. 5 
Fig. 6 
Fig. 7 
Fig. 8 

Fig. 9 

Fig. 10 
Fig. 11 
Fig. 12 
Fig. 13 
Fig. 14 
Fig. 15 
Fig. 16 
Fig. 17 
18 


Imaginary part of relative permeability versus 
frequency for 0.25 dB magnetic foam. 

Relative magnetic susceptibility versus time for 0.25 
dB magnetic foam. 

Real part of relative permeability versus frequency for 
60 dB magnetic foam. 

Imaginary part of relative permeability versus 
frequency for 60 dB magnetic foam. 

Relative magnetic susceptibility versus time for 60 dB 
magnetic foam. 

Real part of relative permeability versus frequency for 
Nickel Ferrite with /i s = 100, - 1 , r 0 = 20 ns. 

Imaginary part of relative permeability versus 
frequency for Nickel Ferrite with /i s - 100, = 1, r Q 

= 20 ns. 

Relative magnetic susceptibility versus time for Nickel 
Ferrite. 

Co-polarized far zone scattered field versus time for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
0.25 dB dielectric foam sphere with 20 cm radius. 

Co-polarized far zone scattered field versus time for 
0.25 dB magnetic foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
0.25 dB magnetic foam sphere with 20 cm radius. 

Co-polarized far zone scattered field versus time for 
60 dB dielectric foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
60 dB dielectric foam sphere with 20 cm radius. 

Co-polarized far zone scattered field versus time for 
60 dB magnetic foam sphere with 20 cm radius. 

Co-polarized Radar Cross Section versus frequency for 
60 dB magnetic foam sphere with 20 cm radius. 


Fig. 



Fig. 

Fig. 
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19 Co-polarized far zone scattered field versus time for 
Nickel Ferrite sphere with 6.67 m radius using both 
dispersive FDTD and non-dispersive FDTD. 

20 Co-polarized Radar Cross Section versus frequency for 
Nickel Ferrite sphere with 6.67 m radius using both 
dispersive FDTD and non-dispersive FDTD. 




Real part of 

0.25 d 


o 

CO 



{ J 7V } ay 


equency (GHz) 



Imaginary pai 

0.25 



o 

lO 

05 


co 

CO 

r- 

C\2 

co 

o 

o 

O 

— - 


c \? 

CV2 

CO 

CO 

o 

d 

d 

o 

d 

d 

o 

o 

d 


{ *r1 | inj 


0.45 




agnetic susce 



Time (ns 









o 

CO 


cs o 







o o o 

O 

o 

o 

o 

O 

O 

o 

o° 

CD C\2 

CO 


o 

CD 


OD 


o 

1 ^ 

1 

1 

cu 

1 

{ 

CO CO 
1 1 

J ri } mj 

1 

1 

lO 

1 

CD 

1 


ency (GHz) 







Magnetic susceptibility, y m , versus time 



i 


(eOI X) mX 


Time (ns 







Magnetic susceptibility, y IT1 , versus time 

Nickel Ferrite 









Co-polarized far zone scattered field 
0.25 dB dielectric foam sphere, 20 cm radius 



Time (ns 










02 

tO 

r- 

o 

cv 

tO 

CO 

o 

CO 

ID 

CO 

CO 

CO 

o 

co 

ID 

09 

o 

CO 

to 

00 

— 

d 

-4 

— I 

d 

o 

d 

c 

d 

o 

o 

d 


(m/^) P[9IJ 0IJ^09[9 9UOZ JlSj 





ection 



Frequency (GHz) 










Time (ns 


















( 2 ui/ap) uonoss 


GO OO 







Co-polarized far zone scattered field 

Nickel Ferrite sphere, 20 cm radius 


'O 



Time 








arized Radar Cross Section 
Ferrite sphere, 6.67 m radius 



equency (MHz) 




